#include "IsoformsEnumerator.h"
// this function returns the next start of the low coverage region.
int nextBeach(CJunctions& ocean, int pos, bool forward)
{
    vector<CJunction>::iterator j = ocean.nextJunction(pos, pos);
    int geneStart = 0;
    int geneEnd = std::numeric_limits<int>::max() - 1; //fake number to the end
    if (forward) {
        return (j == ocean.JunctionsV.end() ? geneEnd + 1 : j->start );
    } else {
        return (j == ocean.JunctionsV.end() ? geneStart - 1 : j->end );
    }
}
////// The functions check if there is a bridge (in, out) the rangeStart and rangeEnd //////
// ptr specifys that we are searching bridge connnected in or conneected out
bool hasBridgeInRange(CJunctions& bridges, int rangeStart, int rangeEnd, bool bridgeOut)
{
    int rangeS = min(rangeStart, rangeEnd);
    int rangeE = max(rangeStart, rangeEnd);
    vector<CJunction>* bridgesV = &bridges.JunctionsV;
    for (vector<CJunction>::iterator it = bridgesV->begin();
            it != bridgesV->end(); it++) {
        int bridgePosi;
        if ( bridges.forward == bridgeOut ) {
            bridgePosi = it->start;
        } else {
            bridgePosi = it->end;
        }
        if (rangeS <= bridgePosi && bridgePosi <= rangeE ) {
            return(true);
        }
    }
    return(false);
}
bool hasBridgeOut(CJunctions& bridges, int rangeStart, int rangeEnd)
{
    return(hasBridgeInRange(bridges, rangeStart, rangeEnd, true));
}

bool hasBridgeIn(CJunctions& bridges, int rangeStart, int rangeEnd)
{
    return(hasBridgeInRange(bridges, rangeStart, rangeEnd, false));
}

// This private function checks if any junction fall in front of a given possition.
// This is called by getPossibleIsoformStarts
bool hasBridgeConntect2(CJunctions& bridges, int loci, int detectRange)
{
    int BEACH_ZONE = bridges.forward ? detectRange : -detectRange;
    // if there is a bridge in the range, return true
    return(hasBridgeIn(bridges, loci, loci + BEACH_ZONE));
}

//////////////////////////////////////////////////////////////////////////
CIsoformsEnumerator::CIsoformsEnumerator(int geneStart, int geneEnd,\
        bool forward, CJunctions* js, CJunctions* ocean)
{
    this->bridges = js;
    this->ocean = ocean;
    this->geneStart = min(geneStart, geneEnd);
    this->geneEnd = max(geneStart, geneEnd);
    this->forward = forward;
    this->clear();
    this->sortJunctions_oceans();
}

CIsoformsEnumerator::~CIsoformsEnumerator(void)
{
}

void CIsoformsEnumerator::clear(void)
{
    this->isoform.exons.clear();
    this->isoforms.clear();
    this->IsoformsKey.clear();
}

void CIsoformsEnumerator::sortJunctions_oceans(void)
{
    this->ocean->sortJunctions(this->forward);
    this->bridges->sortJunctions(this->forward);
}

int CIsoformsEnumerator::enumerateIsoforms(int maxExons)
{
    if (this->forward) {
        return(enumerateIsoforms(this->geneStart, maxExons, this->geneEnd - this->geneStart + 1));
    } else {
        return(enumerateIsoforms(this->geneEnd, maxExons, this->geneEnd - this->geneStart + 1));
    }
}

// The core recursive function to enumerate isoforms
int CIsoformsEnumerator::enumerateIsoforms(int pos, int maxExons, int maxLength)
{
    if (maxExons <= 0 || maxLength <= 0 || pos < geneStart || pos > geneEnd ) {
        if (this->isoform.exons.size() > 0) {
            string key = isoformTag(this->isoform);
            if (this->IsoformsKey.find(key) == this->IsoformsKey.end()) {
                this->IsoformsKey.insert(key);
                this->isoforms.push_back(this->isoform);
            } else {
                printf("\rDuplicate%s", key.c_str());
            }
        }
        return(0);
    }

    // TODO: check if a bridge goes to the ocean.
    // get next zero expressed region (according to direction)
    int nextbeach = this->nextBeach(pos);
    // add the last exon to the no-coverage region, only if there are no junction (splicing) right away.
    if (this->forward) {
        // Try next junction
        for (vector<CJunction>::iterator j = this->bridges->nextJunction(pos, pos);
                j != this->bridges->JunctionsV.end(); j++) {
            int exonBondary = (j->start < this->geneEnd ? j->start : this->geneEnd);
            if (nextbeach < exonBondary) {
                break;
            }
            pair<int, int> exon(pos, exonBondary);
            if (this->isoform.addExon(exon, this->forward)) {
                this->enumerateIsoforms(j->end, maxExons - 1, maxLength - (j->start - pos + 1));
                this->isoform.deleteExon();
            }
        }
        // if this exon has any junction out, don't end in this exon.
        if (!hasBridgeOut(*(this->bridges), pos, nextbeach)) {
            addLastExon(pos, nextbeach - 1);
        }
    } else {
        // Try next junction
        for (vector<CJunction>::iterator j = this->bridges->nextJunction(pos, pos);
                j != this->bridges->JunctionsV.end(); j++) {
            int exonBondary = (j->end > this->geneStart ? j->end : this->geneStart);
            if (nextbeach > exonBondary) {
                break;
            }
            pair<int, int> exon(exonBondary, pos);
            if (this->isoform.addExon(exon, this->forward)) {
                this->enumerateIsoforms(j->start, maxExons - 1, maxLength - (pos - j->end + 1));
                this->isoform.deleteExon();
            }
        }
        // if this exon has any junction out, don't end in this exon.
        if (!hasBridgeOut(*(this->bridges), pos, nextbeach)) {
            addLastExon(nextbeach + 1, pos);
        }
    }
    return(0);
}

int CIsoformsEnumerator::addLastExon(int start, int end)
{
    /* if(this->isoform.exons.size() == 0) {
    	return(0); // no single exon in the middle
    } */

    pair<int, int> exon(start, end);
    if (this->isoform.addExon(exon, this->forward)) {
        string key = isoformTag(this->isoform);
        if (this->IsoformsKey.find(key) == this->IsoformsKey.end()) {
            this->IsoformsKey.insert(key);
            this->isoforms.push_back(this->isoform);
        } else {
            // printf("duplicate"); // remove duplicated isoforms
        }
        this->isoform.deleteExon();
    }
    return(0);
}

int CIsoformsEnumerator::nextBeach(int pos)
{
    vector<CJunction>::iterator j = this->ocean->nextJunction(pos, pos);
    if (forward) {
        return (j == this->ocean->JunctionsV.end() ? this->geneEnd + 1 : j->start );
    } else {
        return (j == this->ocean->JunctionsV.end() ? this->geneStart - 1 : j->end );
    }
}


/* This function find the possible isoform starts by looking for
 * the end of ocean that has no bridge around.
 * Because the region is starting to have coverage,
 * it may be a sign for the start of the transcription. */
vector<int> getPossibleIsoformStarts(int initialStart, CJunctions& bridges, CJunctions& oceans)
{
    vector<int> possibleIstarts;
    possibleIstarts.push_back(initialStart);
    // check through every ocean
    for (vector<CJunction>::iterator it = oceans.JunctionsV.begin();
            it != oceans.JunctionsV.end(); it++) {
        int checkLoci = (bridges.forward ? it->end: it->start);
        bool bValidStart = bridges.forward ? \
                           (initialStart <= checkLoci) :(initialStart >= checkLoci);

        int BEACH_ZONE = 100; // DEFAULT VALUE
        //int BEACH_ZONE = MIN_MRNA_SIZE; // DEFAULT VALUE
        if (MERGE_ALL_ISOFORMS_WITH_THE_SAME_INNTER_EXONS) {
            BEACH_ZONE = abs(nextBeach(oceans, checkLoci, bridges.forward) - checkLoci);
        }
        bValidStart &= !hasBridgeConntect2(bridges, checkLoci, BEACH_ZONE);
        if (bValidStart) {
            int possibleInsoformStart = (bridges.forward ? checkLoci + 1 : checkLoci - 1);
            possibleIstarts.push_back(possibleInsoformStart );
        }
    }
    LOG_INFO("Info %d: There are %d isoform start positions\n",\
             FINE_LOG, (int)possibleIstarts.size());
    return(possibleIstarts);
}

// Get all coding seqeunces from many isoforms
set<string> getCDSs(vector<CTIsoform>& isoforms, const char* contig)
{
    set<string> CDSs;
    int i = 0;
    for (vector<CTIsoform>::iterator it = isoforms.begin();
            it != isoforms.end(); it++, i++) {
        char mRNABuf[MAX_MRNA_SIZE];
        getMRNASeq(*it, contig, mRNABuf);
        getCDS(mRNABuf, CDSs);
    }
    return(CDSs);
}

string isoformTag(CTIsoform& isoform)
{
    string tag;
    char buffer[MAX_LINE];
    unsigned int i;
    for (i = 0; i < isoform.exons.size(); i++) {
        sprintf(buffer, ".%d_%d", isoform.exons.at(i).first, isoform.exons.at(i).second);
        tag = tag + string(buffer);
    }
    return(tag);
}

CIsoformInfo getIsoformInfo(string geneName, const char* isoformId,\
                            const char* contig, CTIsoform& isoform)
{
    int mRNALength = 0;
    for ( int i = 0; i < (int)isoform.exons.size(); i++) {
        mRNALength += abs(isoform.exons.at(i).second - isoform.exons.at(i).first + 1);
    }
    char* mRNABuf = new char[mRNALength + 1];
    mRNALength = getMRNASeq(isoform, contig, mRNABuf);
    CIsoformInfo isoformInfo(isoform);
    isoformInfo.isoformId = isoformId;
    isoformInfo.seq = string(mRNABuf);
    isoformInfo.frame = getTheLongestORF(mRNABuf);
    isoformInfo.CDSscore = evaluateFrame(isoformInfo.frame, isoform.exons);
    isoformInfo.aaSize = (isoformInfo.frame.second - isoformInfo.frame.first)/3;
    isoformInfo.geneName = geneName;
    delete [] mRNABuf;
    return(isoformInfo);
}

// return the isoform has passed the evaluation and print or not
bool filterIsoform(CIsoformInfo& isoformInfo)
{
    // This is where to change the filtering criteria.
    const int MIN_NO_OF_EXONS = 3;
    const int MIN_LENGTH_OF_ISOFORM = 100;
    return((isoformInfo.aaSize >= MIN_LENGTH_OF_ISOFORM
            && (int)isoformInfo.exons.size() >= MIN_NO_OF_EXONS)
           || isoformInfo.CDSscore >= MIN_FRAME_SCORE);
}

//////The following functions are for test purpose//////
// reverse the junction and beach direction (for test purpose only)
void reverseJunctionsAndBeachDirection(CIsoformsEnumerator &e)
{
    e.forward = !e.forward;
    e.bridges->forward = e.forward;
    e.bridges->sortJunctions(e.forward);
    e.ocean->forward = e.forward;
    e.ocean->sortJunctions(e.forward);
}

// This function set up simple junctions and oceans to see
// the forward and backward isoform enumeation is as expectations.
bool testCIsoformsEnumerator(void)
{
    // Setup
    bool pass = true;
    //const char* contig = "012345678901234567890123456789";
    const char* contig = "ACGTACGTACGTACGTACGTACGTACGTACGTACGT";
    const int MAX_EXON = 3;
    int startPos = 0;
    int endPos = (int)strlen(contig) - 1;
    bool forward = true;
    CJunctions js(string("junction"), forward);
    js.addJunction(CJunction(10,50));
    js.addJunction(CJunction(10,100));
    js.addJunction(CJunction(200,300));
    js.sortJunctions(forward);

    CJunctions ocean(string("ocean"), forward);
    ocean.addJunction(CJunction(70,80));
    ocean.addJunction(CJunction(210,290));
    ocean.sortJunctions(forward);

    CIsoformsEnumerator e(startPos, endPos, forward, &js, &ocean);
    // test forward
    e.enumerateIsoforms(startPos, MAX_EXON, (endPos - startPos + 1));
    int EXP_NO_ISOFORMS = 2;
    pass = assert_int_equal(EXP_NO_ISOFORMS, e.isoforms.size());
    printIsoformsExons(cout, e.isoforms);
    // printIsoforms(cout, e.isoforms);
    // test backward
    e.clear();
    reverseJunctionsAndBeachDirection(e); // for convenience, reverse the junctions and ocean.
    e.enumerateIsoforms(endPos, MAX_EXON, (endPos - startPos + 1));
    EXP_NO_ISOFORMS = 1;
    pass = assert_int_equal(EXP_NO_ISOFORMS, e.isoforms.size());
    printIsoformsExons(cout, e.isoforms);
    // printIsoforms(cout, e.isoforms);
    return(pass);
}

bool testPossibleStartSites(void)
{
    // The test is for MERGE_ALL_ISOFORMS_WITH_THE_SAME_INNTER_EXONS = true
    bool pass = true;
    bool forward = true;
    int startPos = 0, endPos = 2000;
    CJunctions js(string("junction"), forward);
    js.addJunction(CJunction(200,400));
    js.addJunction(CJunction(700,800));
    js.addJunction(CJunction(1000,1600));
    js.sortJunctions(forward);

    CJunctions ocean(string("ocean"), forward);
    ocean.addJunction(CJunction(210,390));
    ocean.addJunction(CJunction(701,790));
    ocean.addJunction(CJunction(1200,1400));
    ocean.addJunction(CJunction(1700,1800));
    ocean.sortJunctions(forward);
    // Test forward direction
    int EXP_START_NO;
    int TRUE_START_NO;
    CIsoformsEnumerator e(startPos, endPos, forward, &js, &ocean);
    vector<int> possibleStarts = getPossibleIsoformStarts(startPos, js, ocean);
    EXP_START_NO = 3;
    TRUE_START_NO = possibleStarts.size();
    printIntVector(possibleStarts);
    pass = pass && assert_int_equal(EXP_START_NO, TRUE_START_NO);;
    // Test reverse direction
    reverseJunctionsAndBeachDirection(e);
    possibleStarts = getPossibleIsoformStarts(endPos, js, ocean);
    EXP_START_NO = 3;
    TRUE_START_NO = possibleStarts.size();
    printIntVector(possibleStarts);
    pass = pass && assert_int_equal(EXP_START_NO, TRUE_START_NO);
    return(pass);
}
